Introduction

This file explains how we ran the statisticla tests on the aperture size measurements. It is based on the Aps dataframe, which we created in the Basic Shell Information.R file. It is worth to redo these steps because we won’t summarise it as much this time. First, we need to source the ‘Farasan Database.R’ file through:

source('Farasan Database.R')

Then we extract the aperture measurements from the Apertures list:

Aps <- data.frame()
f <- data.frame()

for (i in 1:length(Apertures))  {
  for (j in 1:length(Apertures[[i]]))  {
    if ("Strombus.fasciatus..............Born.1778..Aperture" %in% names(Apertures[[i]][[j]]) == FALSE) {
      next
    }
    f <-
      data.frame(c(as.character(names(Apertures[i]))), c(as.character(names(Apertures[[i]][j]))), Apertures[[i]][[j]][["Strombus.fasciatus..............Born.1778..Aperture"]])
    colnames(f) <- c('Site', 'Sample', 'Aperture')
    Aps <- rbind(Aps, f)
  }
}
rm(f, i, j, Areas, Sites, Species)

To get a simple list of each measurement grouped by site and bulk sample number:

Aps[1:3,1:3]
##     Site Sample Aperture
## 1 KM1048    926    24.02
## 2 KM1048    926    18.13
## 3 KM1048    926    20.07

General aperture statistics by areas

To look at the general statistics for each areas, we added a column to Aps containing the respective bay areas and summarise the dataframe using the summarySE function:

Aps$Area <- c(rep("Khur Maadi",4122), rep("Janaba West",6171),rep("Janaba East",5169))

source("summarySE.R")

Aps_area <-summarySE(na.omit(Aps), measurevar = "Aperture", groupvars = c("Area"))
##          Area    N Aperture_mean Aperture_median       sd
## 1 Janaba East 5169      19.67116           19.60 1.483827
## 2 Janaba West 6171      21.39223           21.30 1.648272
## 3  Khur Maadi 4122      22.87554           22.85 1.821071

One-way anova

Doing an one-way analysis of variance (anova) will show whether there are statistically significant differences between the means of each bay area.

Aps_area_aov <- aov(Aperture ~ Area, data = Aps)
summary(Aps_area_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)    
## Area            2  23880   11940    4415 <2e-16 ***
## Residuals   15459  41808       3                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post-hoc tukey test

We then applied a post-hoc tukey test, to understand between which areas the statistically significant differences lie.

TukeyHSD(Aps_area_aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Aperture ~ Area, data = Aps)
## 
## $Area
##                             diff      lwr      upr p adj
## Janaba West-Janaba East 1.721064 1.648385 1.793743     0
## Khur Maadi-Janaba East  3.204381 3.123888 3.284873     0
## Khur Maadi-Janaba West  1.483317 1.405778 1.560856     0

Graph of size distributions by area (Fig. 6a)

To illustrate these differences, our paper used raincloud plots, as well as other design-related packages, such as viridis and cowplot.

library(ggplot2)
library(cowplot)
library(viridis)
source('R_rainclouds.R')

Fig6a <- ggplot(Aps, aes(
  x = Area,
  y = Aperture,
  fill = Area,
  colour = Area
)) +
  geom_flat_violin(
    position = position_nudge(x = +0.155, y = 0),
    adjust = 0.75,
    trim = FALSE,
    alpha = .99
  ) +
  geom_point(position = position_jitter(width = .1),
             size = 0.05,alpha=0.2
             
  ) +
  geom_boxplot(
    position = position_nudge(x = +0.155, y = 0),
    aes(x = Area, y = Aperture,colour=Area),
    outlier.shape = NA,
    alpha = 0.3,
    width = .1,
    colour = "BLACK",show.legend = FALSE
  ) +
  ylab('Aperture size') + xlab('') + coord_flip() + guides(fill = FALSE, colour = FALSE) +
  scale_color_manual(values=c(viridis_pal()(60)[10],viridis_pal()(6)[4],viridis_pal()(60)[59]))+
  scale_fill_manual(values=c(viridis_pal()(60)[10],viridis_pal()(6)[4],viridis_pal()(60)[59]))+
  theme_cowplot()+
  theme(axis.text.y.left = element_text(vjust=-2))+theme(axis.ticks.y=element_blank())
Fig6a

General aperture statistics by sites

To do the same tests between sites, we have to take similar steps, only this time we are also using the Site column to group aperture measurements.

##Kruskal Wallis The same data on site level reveals that there are variations in sample size, normal distributions and variance between the sites, which are all assumptions for the ANOVA test used above. We will thus switch to the non-parametric Wilcoxon rank-sum test and use the pairwise.wilcox.test() function

PW_Aps <- pairwise.wilcox.test(Aps$Aperture, Aps$Site,  p.adjust.method = "bonferroni")
PW_Aps
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  Aps$Aperture and Aps$Site 
## 
##          KM1048  KM1050  KM1051  KM1052  KM1053  KM1054  KM1056  KM1057 
## KM1050   1.00000 -       -       -       -       -       -       -      
## KM1051   1.00000 0.60812 -       -       -       -       -       -      
## KM1052   1.00000 0.02953 1.00000 -       -       -       -       -      
## KM1053   1.00000 1.00000 1.00000 1.00000 -       -       -       -      
## KM1054   1.00000 1.00000 1.00000 0.60836 1.00000 -       -       -      
## KM1056   1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 -       -      
## KM1057   1.00000 1.00000 1.00000 0.15740 1.00000 1.00000 1.00000 -      
## KM1304   1.00000 0.01776 1.00000 1.00000 1.00000 1.00000 1.00000 0.10905
## KM1307   1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## KM1313   1.00000 1.00000 0.21637 0.00616 1.00000 1.00000 1.00000 1.00000
## KM1317   1.00000 1.00000 1.00000 0.67309 1.00000 1.00000 1.00000 1.00000
## KM1328   1.00000 1.00000 1.00000 0.19374 1.00000 1.00000 1.00000 1.00000
## KM1330   1.00000 1.00000 1.00000 0.11719 1.00000 1.00000 1.00000 1.00000
## KM1335   1.00000 1.00000 1.00000 0.20132 1.00000 1.00000 1.00000 1.00000
## KM1336   1.00000 1.00000 1.00000 0.48361 1.00000 1.00000 1.00000 1.00000
## JW1705   1.00000 7.7e-06 1.00000 1.00000 1.00000 0.52703 1.00000 0.00011
## JW1727_A 1.00000 < 2e-16 0.00331 1.00000 3.7e-16 0.00090 1.00000 < 2e-16
## JW1727_B 1.00000 < 2e-16 0.00051 1.00000 < 2e-16 0.00047 1.00000 < 2e-16
## JW1807   1.00000 1.2e-09 1.00000 1.00000 1.00000 1.00000 1.00000 < 2e-16
## JW1864   1.00000 < 2e-16 1.5e-05 0.57737 < 2e-16 0.00013 1.00000 < 2e-16
## JW2298   1.00000 < 2e-16 0.04956 1.00000 9.3e-12 0.00231 1.00000 < 2e-16
## JW3120   1.00000 < 2e-16 1.00000 1.00000 2.0e-07 0.02067 1.00000 < 2e-16
## JE0003   1.00000 0.00153 0.19107 1.00000 0.00407 0.19211 1.00000 0.00202
## JE0004   1.00000 < 2e-16 < 2e-16 1.2e-15 < 2e-16 3.4e-10 1.00000 < 2e-16
## JE0078   1.00000 < 2e-16 < 2e-16 4.1e-10 < 2e-16 2.5e-08 1.00000 < 2e-16
## JE0087   1.00000 < 2e-16 5.0e-16 1.1e-06 < 2e-16 3.9e-07 1.00000 < 2e-16
## JE5642   1.00000 < 2e-16 5.3e-08 0.00568 < 2e-16 2.1e-05 1.00000 < 2e-16
##          KM1304  KM1307  KM1313  KM1317  KM1328  KM1330  KM1335  KM1336 
## KM1050   -       -       -       -       -       -       -       -      
## KM1051   -       -       -       -       -       -       -       -      
## KM1052   -       -       -       -       -       -       -       -      
## KM1053   -       -       -       -       -       -       -       -      
## KM1054   -       -       -       -       -       -       -       -      
## KM1056   -       -       -       -       -       -       -       -      
## KM1057   -       -       -       -       -       -       -       -      
## KM1304   -       -       -       -       -       -       -       -      
## KM1307   1.00000 -       -       -       -       -       -       -      
## KM1313   0.00520 1.00000 -       -       -       -       -       -      
## KM1317   0.58399 1.00000 1.00000 -       -       -       -       -      
## KM1328   0.11911 1.00000 1.00000 1.00000 -       -       -       -      
## KM1330   0.10090 1.00000 1.00000 1.00000 1.00000 -       -       -      
## KM1335   0.20101 1.00000 1.00000 1.00000 1.00000 1.00000 -       -      
## KM1336   0.47794 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 -      
## JW1705   1.00000 1.00000 6.4e-06 1.00000 0.01718 0.06750 0.13757 0.67647
## JW1727_A 1.00000 0.00028 < 2e-16 1.1e-10 8.1e-12 2.8e-10 2.4e-13 1.9e-11
## JW1727_B 1.00000 9.5e-05 < 2e-16 1.4e-11 3.0e-12 7.0e-11 7.1e-14 2.3e-12
## JW1807   1.00000 1.00000 8.2e-08 1.00000 0.01638 0.10223 0.07202 0.84577
## JW1864   1.00000 5.6e-06 < 2e-16 1.9e-14 5.5e-15 4.5e-13 < 2e-16 2.4e-15
## JW2298   1.00000 0.00087 < 2e-16 5.0e-09 2.7e-09 4.1e-09 1.2e-10 1.4e-09
## JW3120   1.00000 0.04652 < 2e-16 7.9e-06 3.0e-07 1.8e-06 1.1e-07 2.5e-06
## JE0003   1.00000 0.04016 0.00086 0.00382 0.02436 0.00269 0.00533 0.00465
## JE0004   3.6e-16 3.3e-15 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16
## JE0078   9.2e-10 1.2e-12 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16
## JE0087   5.3e-06 6.1e-11 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16
## JE5642   0.01192 4.5e-07 < 2e-16 6.5e-15 4.1e-15 3.6e-13 < 2e-16 1.0e-15
##          JW1705  JW1727_A JW1727_B JW1807  JW1864  JW2298  JW3120  JE0003 
## KM1050   -       -        -        -       -       -       -       -      
## KM1051   -       -        -        -       -       -       -       -      
## KM1052   -       -        -        -       -       -       -       -      
## KM1053   -       -        -        -       -       -       -       -      
## KM1054   -       -        -        -       -       -       -       -      
## KM1056   -       -        -        -       -       -       -       -      
## KM1057   -       -        -        -       -       -       -       -      
## KM1304   -       -        -        -       -       -       -       -      
## KM1307   -       -        -        -       -       -       -       -      
## KM1313   -       -        -        -       -       -       -       -      
## KM1317   -       -        -        -       -       -       -       -      
## KM1328   -       -        -        -       -       -       -       -      
## KM1330   -       -        -        -       -       -       -       -      
## KM1335   -       -        -        -       -       -       -       -      
## KM1336   -       -        -        -       -       -       -       -      
## JW1705   -       -        -        -       -       -       -       -      
## JW1727_A < 2e-16 -        -        -       -       -       -       -      
## JW1727_B < 2e-16 1.00000  -        -       -       -       -       -      
## JW1807   1.00000 < 2e-16  < 2e-16  -       -       -       -       -      
## JW1864   < 2e-16 0.00316  1.00000  < 2e-16 -       -       -       -      
## JW2298   3.7e-10 1.00000  1.00000  < 2e-16 0.35921 -       -       -      
## JW3120   9.5e-05 1.4e-05  1.9e-07  2.0e-12 < 2e-16 0.58862 -       -      
## JE0003   0.01027 1.00000  1.00000  0.01014 1.00000 1.00000 0.38053 -      
## JE0004   < 2e-16 < 2e-16  < 2e-16  < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.63341
## JE0078   < 2e-16 < 2e-16  < 2e-16  < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.00000
## JE0087   < 2e-16 < 2e-16  < 2e-16  < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.00000
## JE5642   < 2e-16 8.3e-08  0.00032  < 2e-16 0.00197 3.8e-06 1.5e-15 1.00000
##          JE0004  JE0078  JE0087 
## KM1050   -       -       -      
## KM1051   -       -       -      
## KM1052   -       -       -      
## KM1053   -       -       -      
## KM1054   -       -       -      
## KM1056   -       -       -      
## KM1057   -       -       -      
## KM1304   -       -       -      
## KM1307   -       -       -      
## KM1313   -       -       -      
## KM1317   -       -       -      
## KM1328   -       -       -      
## KM1330   -       -       -      
## KM1335   -       -       -      
## KM1336   -       -       -      
## JW1705   -       -       -      
## JW1727_A -       -       -      
## JW1727_B -       -       -      
## JW1807   -       -       -      
## JW1864   -       -       -      
## JW2298   -       -       -      
## JW3120   -       -       -      
## JE0003   -       -       -      
## JE0004   -       -       -      
## JE0078   < 2e-16 -       -      
## JE0087   < 2e-16 < 2e-16 -      
## JE5642   < 2e-16 4.5e-08 1.00000
## 
## P value adjustment method: bonferroni

Graph of size distributions by site (Fig. 6b)

Similar to the graph above, we illustrate the site differences using raincloud plots, as well as other design-related packages, such as viridis and cowplot. This time, the sites are coloured by bay area, however they are in the order of largst shells (top) to smallest shells (bottom). The general trend, however, remains and sites usually stick to their bay area.

library(ggplot2)
library(cowplot)
library(viridis)
source('R_rainclouds.R')

 Fig6b <- ggplot(Aps, aes(
    x =   reorder(Site, Aperture, FUN = mean,order=TRUE),
    y = Aperture,
    fill = Area,
    colour = Area, order=Aperture
  )) +
    geom_flat_violin(
      position = position_nudge(x = +0.2, y = 0),
      adjust = 0.75,
      trim = FALSE,
      alpha = .99
    ) +
    geom_point(position = position_jitter(width = 0.1),
               size = 0.5,alpha=0.3
    ) +
    geom_boxplot(
      position = position_nudge(x = +0.2, y = 0),
      aes(x = Site, y = Aperture,colour=Area),
      outlier.shape = NA,
      alpha = 0.3,
      width = .1,
      colour = "BLACK"
    ) +
    ylab('Aperture size') + xlab("") + coord_flip() +guides(fill = guide_legend(reverse = TRUE),color = guide_legend(reverse = TRUE))+
    scale_color_manual(values=c(viridis_pal()(3)[1],viridis_pal()(6)[4],viridis_pal()(60)[59]))+
    scale_fill_manual(values=c(viridis_pal()(3)[1],viridis_pal()(6)[4],viridis_pal()(60)[59]))+    theme_cowplot()+ ggtitle("")+
    theme(axis.text.y.left = element_text(vjust=-0))+
    theme(axis.ticks.y=element_blank())+ 
    theme(legend.position="right")
  
  Fig6b

Finally, we want to integrate both raincloud plots (Areas and Sites) into one figure. For this we use the patchwork package.

library(patchwork)
Fig6 <- Fig6a+Fig6b

Fig6 + plot_annotation(
  title = 'Aperture sizes',
  subtitle = 'Differences between bay area on the left and differences between sites on the right')

setwd("~/Farasan-Shell-Data-master/Figures and supplementary material")#Remember to change the file path
ggsave("Figure_6.png", height = 200, width = 300, units = "mm")